1 Background

tmb <- readRDS("depends/tmb.rds")
aghq <- readRDS("depends/aghq.rds")
tmbstan <- readRDS("depends/tmbstan.rds")

All of the possible parameter names are as follows:

unique(names(tmb$fit$obj$env$par))
##  [1] "beta_rho"             "beta_alpha"           "beta_lambda"          "beta_anc_rho"         "beta_anc_alpha"       "logit_phi_rho_x"      "log_sigma_rho_x"      "logit_phi_rho_xs"     "log_sigma_rho_xs"     "logit_phi_rho_a"     
## [11] "log_sigma_rho_a"      "logit_phi_rho_as"     "log_sigma_rho_as"     "log_sigma_rho_xa"     "u_rho_x"              "us_rho_x"             "u_rho_xs"             "us_rho_xs"            "u_rho_a"              "u_rho_as"            
## [21] "logit_phi_alpha_x"    "log_sigma_alpha_x"    "logit_phi_alpha_xs"   "log_sigma_alpha_xs"   "logit_phi_alpha_a"    "log_sigma_alpha_a"    "logit_phi_alpha_as"   "log_sigma_alpha_as"   "log_sigma_alpha_xa"   "u_alpha_x"           
## [31] "us_alpha_x"           "u_alpha_xs"           "us_alpha_xs"          "u_alpha_a"            "u_alpha_as"           "u_alpha_xa"           "OmegaT_raw"           "log_betaT"            "log_sigma_lambda_x"   "ui_lambda_x"         
## [41] "log_sigma_ancrho_x"   "log_sigma_ancalpha_x" "ui_anc_rho_x"         "ui_anc_alpha_x"       "log_sigma_or_gamma"   "log_or_gamma"

2 Time taken

data.frame(
  "TMB" = tmb$time,
  "aghq" = aghq$time,
  "tmbstan" = tmbstan$time
)

3 Histogram

histogram_plot("beta_anc_rho")

4 KS

4.1 Individual parameters

ks_helper <- function(par) to_ks_df(par) %>% ks_plot(par = par)

ks_helper("beta")

ks_helper("logit")

ks_helper("log_sigma")

ks_helper("u_rho_x")

ks_helper("u_rho_xs")

ks_helper("us_rho_x")

ks_helper("us_rho_xs")

ks_helper("u_rho_a")

ks_helper("u_rho_as")

ks_helper("u_alpha_x")

ks_helper("u_alpha_xs")

ks_helper("us_alpha_x")

ks_helper("us_alpha_xs")

ks_helper("u_alpha_a")

ks_helper("u_alpha_as")

ks_helper("u_alpha_xa")

ks_helper("ui_anc_rho_x")

ks_helper("ui_anc_alpha_x")

ks_helper("log_or_gamma")

4.2 Summary table

ks_summary <- lapply(unique(names(tmb$fit$obj$env$par)), function(x) {
  to_ks_df(x) %>%
    group_by(method) %>%
    summarise(ks = mean(ks), par = x)
}) %>%
  bind_rows() %>%
  pivot_wider(names_from = "method", values_from = "ks") %>%
  rename(
    "Parameter" = "par",
    "KS(aghq, tmbstan)" = "aghq",
    "KS(TMB, tmbstan)" = "TMB",
  )

ks_summary %>%
  gt::gt() %>%
  gt::fmt_number(
    columns = starts_with("KS"),
    decimals = 3
  )
Parameter KS(aghq, tmbstan) KS(TMB, tmbstan)
beta_rho 0.182 0.188
beta_alpha 0.117 0.117
beta_lambda 0.057 0.058
beta_anc_rho 0.111 0.110
beta_anc_alpha 0.097 0.057
logit_phi_rho_x 0.391 0.626
log_sigma_rho_x 0.581 0.656
logit_phi_rho_xs 0.270 0.638
log_sigma_rho_xs 0.839 0.646
logit_phi_rho_a 0.828 0.528
log_sigma_rho_a 0.438 0.557
logit_phi_rho_as 0.824 0.514
log_sigma_rho_as 0.294 0.542
log_sigma_rho_xa 0.438 0.692
u_rho_x 0.116 0.114
us_rho_x 0.103 0.102
u_rho_xs 0.147 0.146
us_rho_xs 0.101 0.101
u_rho_a 0.179 0.180
u_rho_as 0.125 0.118
logit_phi_alpha_x 0.317 0.571
log_sigma_alpha_x 0.709 0.592
logit_phi_alpha_xs 0.279 0.587
log_sigma_alpha_xs 0.686 0.613
logit_phi_alpha_a 0.791 0.533
log_sigma_alpha_a 0.290 0.540
logit_phi_alpha_as 0.749 0.511
log_sigma_alpha_as 0.237 0.540
log_sigma_alpha_xa 0.706 0.653
u_alpha_x 0.091 0.086
us_alpha_x 0.091 0.091
u_alpha_xs 0.078 0.069
us_alpha_xs 0.116 0.114
u_alpha_a 0.126 0.124
u_alpha_as 0.086 0.080
u_alpha_xa 0.073 0.072
OmegaT_raw 0.192 0.548
log_betaT 0.174 0.679
log_sigma_lambda_x 0.556 0.770
ui_lambda_x 0.205 0.205
log_sigma_ancrho_x 0.714 0.524
log_sigma_ancalpha_x 0.706 0.616
ui_anc_rho_x 0.066 0.071
ui_anc_alpha_x 0.132 0.134
log_sigma_or_gamma 0.507 0.521
log_or_gamma 0.075 0.075
 
ggplot(ks_summary, aes(x = `KS(TMB, tmbstan)`, y = `KS(aghq, tmbstan)`)) +
  geom_point(alpha = 0.5) +
  xlim(0, 1) +
  ylim(0, 1) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  labs(x = "KS(aghq, tmbstan)", y = "KS(TMB, tmbstan)") +
  theme_minimal() 

ks_summary %>%
  mutate(`KS difference` = `KS(TMB, tmbstan)` - `KS(aghq, tmbstan)`) %>%
  ggplot(aes(x = `KS difference`)) +
    geom_boxplot(width = 0.5) +
    coord_flip() +
    labs(x = "KS(TMB, tmbstan) - KS(aghq, tmbstan)") +
    theme_minimal()

5 MMD

#' To write!

6 PSIS

#' To write!